An 0(N 2 ) Approximation for Hydrodynamic Interactions in Brownian Dynamics Simulations 
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In the Ermak-McCammon algorithm for Brownian Dynamics, the hydrodynamic interactions (HI) between TV 
spherical particles are described by a 3iV x 3N diffusion tensor. This tensor has to be factorized at each timestep 
with a runtime of 0(N 3 ), making the calculation of the correlated random displacements the bottleneck for 
many-particle simulations. Here we present a faster algorithm for this step, which is based on a truncated 
expansion of the hydrodynamic multi-particle correlations as two-body contributions. The comparison to the 
exact algorithm and to the Chebyshev approximation of Fixman verifies that for bead-spring polymers this 
approximation yields about 95% of the hydrodynamic correlations at an improved runtime scaling of 0(N 2 ) 
and a reduced memory footprint. The approximation is independent of the actual form of the hydrodynamic 
tensor and can be applied to arbitrary particle configurations. This now allows to include HI into large many- 
particle Brownian dynamics simulations, where until now the runtime scaling of the correlated random motion 
was prohibitive. 

PACS numbers: 47.85.Dh, 83.1().Mj 



I. INTRODUCTION 



For the simulation of diffusional processes on the scales of 
polymers or proteins. .Rmwnian Dynamics has proven to be a 
reliable workhorseLLLMJU. This coarse grained method builds 
upon Einstein's microscopic explanation of the random mo- 
tion of colloidal particlesU which had been observed earlier 
by the biologist Robert Brown when studying pollen grains. 
In Einstein's explanation the solvent molecules are replaced 
by a heat bath that models the collisions between the solvent 
molecules and the much larger Brownian particles by cor- 
rectly distributed random forces. Together with the assump- 
tion that the motion of the large observable particles is over- 
damped, i.e., that their velocities relax very fast compared to 
the time steps, used in the simulation, their diffusive motion 
is reproduced!! Replacing the many small solvent molecules 
by a continuum dramatically reduces the computational costs 
compared to an all-atom molecular dynamics simulation, but 
for the price that also all interactions between the Brownian 
particles have to be adapted to include the effects of the now 
continuous solvent. For electrostatic interactions, e.g., the po- 
larizability of the solvent molecules and the redistribution of 
the included ions can be described by an effective shielding 
via the Debye-Hiickel theory. When the solvent molecules 
are not considered explicitly anymore, also so called hydro- 
dynamic interactions (HI) have to be included, which describe 
how the relative motion of the Brownian particles is coupled 
mechanically via the displaced solvent. 

Now there is a fundamental difference between the effect 
of the hydrodynamic coupling onto the direct interactions be- 
tween the Brownian particles and the one onto their random 
motion. An external force acting on one of the particles ac- 
celerates this particle, which in turn leads to a displacement 
of the other particles, too, mediated by the displaced solvent 
molecules. Effectively, with hydrodynamics the effect of the 
external forces is increased and the dynamics is accelerated 
compared to a setup that does not consider hydrodynamics. In 
contrast, the random motion of the Brownian particles models 
the effect of the thermal fluctuations of the solvent molecules, 
the strength of which depends on the temperature of the sys- 



tem. Consequently, for a correct description of the overall 
diffusion the random kicks must have the same strength with 
and without hydrodynamics, i.e., whether the correlation due 
to the displaced solvent is considered or neglected. Conse- 
quently, when hydrodynamic interactions are included in a BD 
simulation, the random displacements still have (up to second 
order) the-same (temperature dependent) magnitudes, but are 
correlated! In a simulation, they now have to be determined 
from a factorization of the diffusion tensor of the complete 
system, which is numerically demanding. 

In a BD simulation_.of N particles using the Ermak- 
McCammon algorithmU, the hydrodynamically correlated 
random displacements are determined via a Cholesky factor- 
ization of the 3N x 3N diffusion tensor, which results in a 
runtime of 0(N 3 ), while all the two-body interactions can be 
calculated in 0(N 2 ) runtime. For a many particle simulation 
with HI, one is therefore spending most of the time evaluating 
the correlated random displacements. For practical applica- 
tions, this limits the number of particles to a few dozens, when 
HI is included, while without HI the dynamics of some thou- 
sand particles could be simulated in the same time. As the 
already approximated direct interactions are most important 
for the dynamics of most biologically interesting association 
and dissociation processes, the time consuming HI is often 
neglectedLJ. For applications to polymers or the dynamics of 
DNA, however, it may be important to explicitly include HI in 
order to reproduce the correct dynamic behaviorlL±jU. 

Considering that HI should on the one hand not be com- 
pletely neglected, but on the other hand is only one out of a 
handful of interactions, there is an urgent need for faster al- 
gorithms to compute the hydrodynamic coupling in the ran- 
dom motion of many particle systems. Fixman proposed to 
use a Chebyshev approximation, which scales as 0(N 2 - 5 ), to 
factorize the diffusion coefficient!!.! A drastically simplified 
approach to HI is the effective mobility model by Heyes et 
al., which changes the diffusion coefficients of the Brownian 
particles according to the local density!—!. Such a model can 
obviously not reproduce the full correlation between the in- 
dividual particles. Based on these ideas, Banchio and Brady 
developed an algorithm for infinite homogenous suspensions 
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that scales as ©(iV 125 log AT) for large Due to its mathe- 
matical complexity, this algorithm only performs well for very 
large numbers of particles. Consequently, when HI is to be 
considered in BD simulations, many current studies either use 
the original Cholesky factorization for smaller systems or Fix- 
man's Chebyshev approximation wittUts better runtime scal- 
ing for larger many particle scenariosLI. 

Here we present a conceptionally different approach to 
tackle the computationally expensive correlation of random 
motion in many-particle simulations. Whereas all previous 
improvements to the Ermak-McCammon algorithm only con- 
sidered the factorization of the hydrodynamic tensor, we ar- 
gue that in typical BD simulations all interactions are approx- 
imative anyhow. Consequently, even with the exact (and very 
time-consuming) HI, the system dynamics will not be exactly 
the same as in an experiment. Thus one may try to approx- 
imate the correlations of the random forces, too, and reduce 
their functional complexity without further perturbing the dy- 
namics in the simulation. The central requirement for such an 
approximation is thus that it may not introduce any systematic 
drift. 

This publication is organized as follows. Before we present 
our truncated expansion ansatz, we shortly review the stan- 
dard Ermak-McCammon algorithm for BD simulations and 
how HI is treated there. We also give a short introduction 
to the Chebyshev expansion of the HI as introduced by Fix- 
man. In section ^l| we compare our ansatz to the well es- 
tablished methods of Ermak and Fixman, respectively. The 
main tests are the dynamics of a bead-spring dimer, which 
was used by Ermak and McCammon to verify their algorithm, 
and the behavior of bead-spring polymers of various lengths. 
The simulation results show that our ansatz reproduces about 
95% of the correlations due to hydrodynamics at a runtime 
which scales quadratic with the bead number. We close with 
a summary and an outlook pointing out potential extensions 
and applications. 



II. EVALUATING HYDRODYNAMIC INTERACTIONS 



covariance, which is described by the corresponding entries of 
D: 



(Ri(At)) = 0, (Ri(At)Rj(M)) = 2D lj At . 



(2) 



The most common forms for the diffusion tensor are-the Os- 
een tensonJ and the Rotne-Prager-Yamakawa tensottilO. For 
these approximations the second term of equation ([I]), which 
describes how the particles are dragged into regions of faster 
diffusion, vanishes. 

For identical spheres of radius a, the 37V x 3N Rotne- 
Prager-Yamakawa hydrodynamic tensor, which couples the 
translational displacements of the beads, consists of the fol- 
lowing 3x3 submatrices Dy , where i and j label two parti- 
cles: 
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This tensor is positive definite for all particle configurations. 
Dickinson et al.U later showed how rotational coupling can 
be included, too. A generalization of the Rotne-Prager- 
Yamakawa tensor for .spheres of different radii was introduced 
by Garcia de la TorreU. 

To determine the random displacements from equation (|2|), 
one needs to find B so that D = BB T . One possible solu- 
tion, which was used by Ermak and McCammon, comes from 
a Cholesky factorization, giving B as an upper (lower) tridiag- 
onal matrix. The random displacements are then determined 
as R = BX from a 37V dimensional vector X of normal dis- 
tributed random numbers. 



A. The Ermak-McCammon algorithm and the Rotne-Prager 
tensor 



B. Fixman's Chebyshev approximation 



In the often mpd Ermak-McCammon algorithm for Brow- 
nian DynamicsUU the total displacement Ar^ of the ith par- 
ticle during a timestep At due to the external forces Fj acting 
on all the particles and the random displacement Ri is given 
by 



Ar,(A<) = V ^lAt 
V ; ^ k B T 



EdDij . 



Ri(At) (1) 



The 37V x 3iV diffusion tensor D = (Dij ) describes the hy- 
drodynamic coupling between the TV particles with their three 
translational degrees of freedom. For the external forces, D 
is used directly, whereas the hydrodynamically correlated ran- 
dom displacements Ri(At) are characterized only indirectly 
by the statistical moments of a vanishing average and a finite 



To avoid the computationally expensive Cholesky fac- 
torization of the hydrodynamic tensor used in the Ermak- 
McCammon algorithm, Fixman suggested in 1986 to approx- 
imate the square root of the diffusion tensor via Chebyshev 
polynomials!.! Fortunately, the explicit calculation of this ma- 
trix is not necessary and the vector of correlated displacements 
can be determined iteratively by a series of matrix-vector mul- 
tiplications up to the order L of the expansion: 



L 

R ~ y^^aixi 

1=0 

x = X 

xi = [d a B + d b I] ■ X 
afj+i = 2[d a T> + d b I] ■ xi - xi- 



(8) 

(9) 
(10) 
(11) 
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The factors d a and db are related to the range of eigenvectors 

[A m in j ^max] Of D as 

and db = (12) 



A 



max f^miri 



A, 



Amaa: A, 



For further details like the evaluation of the expansion coeffi- 
cients ai we refer the readers to, e.g., Press et al.tJ 

To determine the necessary order L of the Chebyshev ap- 
proximation and to control whether the eigenvalue spectrum 
of D is bounded correctly, we followed the procedure of Jen- 
drejack et al£j. They introduced the relative error ef (Ef 
in their notation) derived from the approximated random dis- 
placements R, the uncorrected random numbers X, and the 
diffusion tensor D: 



\R-R-X~DXl 



e / = 



ADA 



(13) 



As R ■ R can be evaluated with negligible additional cost, we 
used ef to monitor the convergence of the Chebyshev iteration 
dl~i|). When Cf remained above the chosen threshold with the 
actual values of L, \ m in< an d A ma x> L was increased by three 
and X m in and X max were recalculated. For the next iteration, 
L = Imax + 3 was used, where l max is the index of the last 
term of the iteration < fFT| ) required to gfti£/ below the chosen 
threshold. Following Jendrejack et al.tJ, ef was set to 10 -3 
if not otherwise noted. For this value the numerical results 
were sufficiently close to the results with the exact Cholesky 
factorization. 



C. The truncated expansion ansatz 

The two methods outlined above focus on factorizing D. 
Compared to the exact Cholesky factorization, the numeri- 
cal approximation of Fixman manages to efficiently get close 
to the exact value of the square root of the diffusion tensor. 
We now present an approximation tailored for practical appli- 
cations of many particle simulations, which algorithmically 
treats the displacements from the external forces and from the 
random forces on an equal footing. 

For such an approximation we start from equation (jlj). Ne- 
glecting the random displacements, the displacement along 
coordinate i during At is 



V ; ^ k B T k B T 1 ' 



(14) 



where we introduced the hydrodynamically corrected effec- 



tive force F, 



eff. 



D, 



(15) 



This reformulation is independent of the actual form of the 
hydrodynamic tensor. 

For the random displacements, we now make an ansatz with 
the same structure, i.e., we also introduce a hydrodynamically 



corrected random force acting on coordinate i which is 
derived from the uncorrelated random forces fj that would act 
on each of the particles in the absence of HI. In this ansatz the 
displacements due to the random forces alone are given as 



Di 



Da At 
k R T 



eff 



(16) 



The structure of this ansatz, which effectively factorizes D 
into individual two-body contributions, is taken from equa- 
tion (P~4|). The scaling factor d takes care that for each indi- 
vidual coordinate its unperturbed diffusion coefficient Da is 
regained! while /Jy allows for different weights when \/D is 
used instead of D in equation (|5j). As we will see later, there 
is an individual scaling factor C{ for each of the coordinates, 
while, due to symmetry reasons, for the coefficients (3ij we 
only need two different values for the diagonal f3u and for the 
off-diagonal ^ with i ^ j, respectively. 

Now the parameters C; and /3y have to be determined such 
that the moments of the correlated displacements (Q) are re- 
produced with only small deviations. For the uncoupled ran- 
dom forces we use 



(/,)=0 and (f lf] ) = ^M_£ S ^ 



(17) 



i.e., they reproduce the mean and covariance of the random 
displacements in the uncorrelated case. 

With correlation, i.e., with hydrodynamic coupling, the 
vanishing average (Ri) is fulfilled straightforwardly from 
(fi) = 0. For the covariance we start from the product of 
Ar.i and Arj according to equation (p~6|): 



An An 



At 



C,C, J2 hi. [,,!),, J ),,]). J) (18) 



k.i 



Inserting (111) leads to the condition 



D ik D jk i 



r irj ) = 2AtCiCj y p ik p jk ~™~ JK = 2DijAt. (19) 



The terms with k ^ I drop from the double sum, because the 
h are uncorrelated. Requiring that the variance of equation 
(g) be reproducedl-l for i = j allows us to determine the nor- 
malization constants C, : 



a 



D U D 



(20) 



ii Ukk 



Without loss of generality, we can set (3a = 1. Then, equation 
dl6| ) reduces to the usual form in the limit of vanishing HI. 
In this case, also Ci = 1 and the displacement Ari (At) can 
be calculated with equation ([IJ) from the sum of the random 
force fi and the external force Fi . 

With [3a = 1, equation ( |20| ) can also be written as 



^ =1 



D 2 

u ik 

DnDkk 



(21) 
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where the term with k = i was taken out of the sum. 



With our ansatz (16) there is no set of coefficients which 



can be determined numerically efficiently, with which equa- 
tions (||) can be fulfilled simultaneously. To determine the 
remaining off-diagonal coefficients ;, we therefore proceed 
by assuming that the hydrodynamic coupling is weak, i.e., that 
the off-diagonal entries D$j of the hydrodynamic tensor are 
much smaller than the individual diffusion coefficients Da. 
Then we can use \/l + e w 1 + e/2 to expand equation (pJJ): 



and Pij according to equation ( |28| ) has the same structure as 
the displacements due to the hydrodynamically coupled exter- 
nal forces in equation (|l4|), resulting in the same runtime scal- 
ing of 0(N 2 ) for both the deterministic and the random con- 
tributions to Ar-j. Consequently, with this algorithm, which 
treats the external and the random forces on equal footing, HI 
can now be included in all those BD simulations where the di- 
rect interparticle forces can be computed — given our approx- 
imation is accurate enough for the chosen application. 



1 

a. 



n 2 

u ik 



2 /_j IK £) £) 

k^i 



(22) 



Thus, the product (CiCj) 1 in equation |l9| can be approxi- 
mated as 



1 



GiCj 



(23) 

i.e., a constant term plus terms that are quadratic and quartic 
in Dij/Du. As we have the same ftu = [3 = 1 for all i, we 
also set all /3y for i ^ j to the same value (3' . Without the 
quartic term, equation ([23]) then becomes 



CiCj 



1 + (N - 1) /3' 2 e 2 



(24) 



with the averaged relative coupling strength e = (Dij/Du). 

Starting directly from equation ([l^) we get the non- 
approximated relation 



1 1 

Ci Cj 



^ PikPjk 



DkkDij 



(25) 



On the rhs of this equation the two terms with k = i and 
k = j can be taken out of the sum. They are independent of 
the hydrodynamic coupling, while all other terms are of first 
order in Dij/Du. With (3=1 this yields: 



1 1 

Ci Cj 



2(3' + (N - 2)(3' 2 e 



(26) 



Comparing equations (Q) and ( [26] ) then gives the quadratic 
equation 

(3' 2 [{N - l)e 2 - (N - 2)e] - 2(3' + 1 = 0, (27) 
which can be solved for (3': 



hi = 



1 - sj\ - [(TV - l)e 2 - (TV — 2)t 
(N — l)e 2 - (N -2)e 



(28) 



We note that according to the above equation is not defined 
in the case of absent HI. However, it converges to = 1/2 
for vanis hing HI, i.e., for e — > 0, the value obtained from 
equation ( P7| ) for e = 0. 

The resulting approximation to the hydrodynamically cou- 
pled random displacements of equation (|l6|) with the normal- 
ization constants C, given by (Eoj) and the two weights (3u = 1 



D. Simulation details 

To evaluate how our ansatz compares to the standard meth- 
ods, we ran BD simulations of bead-spring polymers with 
N = 2. .. 2000 beads of radius a. The diffusion coefficient of 
the individual beads was set to Dq = 1, thus defining the time 
scale. The beads were connected to their direct neighbors in 
the chain by harmonic springs with 



V h (x) = a h (x - L) 2 . 



(29) 



The potential minimum was varied for the dimer simulations 
and set to L = 3a for the polymer simulations. To prevent 
the beads from overlapping, a repulsive harmonic potential 
between all beads. was used analogous to the setup of Ermak 
and McCammonlii]: 



V c (x) = a h (x - 2a) 2 for r < 2a 



(30) 



The spring constant for both interactions was set to the rather 
stiff value of ah = 50 ksT/a 2 , with which we used a conser- 
vatively estimated integration timestep At = 0.001. 

The Cholesky factorization, the eigenvalue calculation, the 
matrix and vector operations as well as the generation of the 
random numbers were performed with the respective subrou- 
tines from the GNU scientific libraryLI. 

For simulations of shorter polymers of N < 70, all simula- 
tion were started with the beads of the polymer aligned along 
the x axis with a mutual separation L = 3. Data analysis 
was started when the polymer had reached its coiled state. 
This transition was observed via the radius of gyration. For 
polymers with N > 70, the equilibration was too slow with 
hydrodynamics included. For every chain length we there- 
fore started a set of simulations without hydrodynamics and a 
longer timestep of At = 0.003 and saved the final positions 
when the polymer had coiled up. These snapshots were then 
used as starting points for simulations with the different forms 
of HI. 



III. COMPARISON OF THE METHODS 

In the following section we evaluate how our truncated ex- 
pansion approximation to HI (TEA-HI) compares to the other 
methods, namely the exact factorization of Ermak and Mc- 
Cammon and the Chebyshev approximation introduced by 
Fixman. For this, we analyzed the correlation coefficient of 
two random displacements in one dimension, the dynamics of 
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(d) 




Dcm 






an 








Geyer 


Ermak 


theo. 


Geyer 


Ermak 


theo 


2 


2.028 


0.7164 


0.7468 


0.7464 


2.43 


2.28 


2.28 


3 


3.003 


0.6530 


0.6665 


0.6665 


2.97 


2.93 


2.93 


4 


4.002 


0.6175 


0.6253 


0.6249 


3.27 


3.27 


3.22 


8 


8.001 


0.5612 


0.5627 


0.5625 


3.59 


3.59 


3.62 


20 


20 


0.5247 


0.5243 


0.5250 


3.85 


3.92 


3.85 


66.7 


66.7 


0.5073 


0.5073 


0.5075 


3.92 


3.91 


3.96 



TABLE I: Comparison of the center-of-mass diffusion coefficient 
Dcm and the rescaled correlation time «d of the orientation of 
a dimer of spherical beads with radius a connected by a spring of 
length L from BD simulations with HI to their respective theoretical 
values. The actual average separation during the simulation is given 
by (d) and was used to calculate the theoretical predictions from ( |32| ) 
and (|33j). The headings "Geyer" and "Ermak" denote the respective 
type of HI that was used in the simulations. 



a dimer, which had been introduced as a test case by Ermak 
and McCammonU, and several static and dynamic properties 
of bead-spring polymers. These comparisons will confirm that 
most of the hydrodynamic interaction is obtained with our ap- 
proximation at a greatly improved runtime scaling. 
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FIG. 1: Average radius of gyration (R%) (diamonds) and end-to- 
end distance (R^e) (points) for polymers of various chain lengths 
TV from simulations where HI was calculated with Fixman's Cheby- 
shev approximation (open symbols) and with our truncated expan- 
sion approximation (filled symbols), respectively. The lines indicate 
the theoretically predicted scaling cx (N — 1) . 



the HI, the analytical values for Dc m and an = Td/td are 



A. Correlation coefficient 

The simplest test is to directly evaluate the correlation co- 
efficient pij of two one-dimensional displacements that are 
hydrodynamically coupled: 



D c 



Pij 



(Ari){A rj 



(31) 



For Dn = D22, D12 = D21 was varied between and 
-Dn and p^j was averaged for correlated displacements using 
both the explicit Cholesky factorization and our approxima- 
tion. For the whole range of D12, our TEA-HI gave values 
for P12, which were by less than 0.1% smaller than with the 
exact factorization. Similar minor deviations were obtained, 
too, for all tested cases with Dn 7^ Z?22- 



B. Dimer dynamics 

In their original workD, Ermak and McCammon verified 
their approach to BD with HI by comparing the numerically 
determined diffusion coefficient of the center of mass of a 
dimer of spherical beads with radius a and separation L, 
Dcm, and the inverse of the rescaled relaxation time of its ori- 
entation, ajj, to analytical results. Only translational coupling 
between the beads was considered with the-Rotne-Prager ten- 
sor. According to Ermak and McCammonU with this form of 
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(32) 
(33) 



Here, d is the distance between the centers of the two beads. 
Dcm is given in units of the diffusion coefficient Do of a 
single bead. The relaxation time td is given in units of To = 
d 2 /D , which is the time that a single bead would need to 
diffuse over the separation d between the two beads of the 
dimer. Without HI, the limiting values are Dcm — 0.5 and 
a D = 4. 

In a spring-bead dimer the average distance between the 
two beads is slightly larger than the length L of the spring 
connecting them. We therefore calculated Dcm and au from 
the observed average separation (d) during the simulation for 
different spring lengths L. As seen in table Q, our results re- 
produce the analytical values for both Dcm and an quite 
well. The hydrodynamic coupling tends to be underestimated 
by less than 7% for touching spheres, where the correlation is 
strongest, and much less for larger separations. Underestimat- 
ing the HI has the consequence that the diffusion coefficient 
Dcm is slightly smaller than the theoretical value and the re- 
laxation time is slightly shorter, resulting in a larger old- For 
comparison, table | also gives the numerical results for Dcm 
and au with the correct factorization of Ermak and McCam- 
mon. Their deviation from the theoretical predictions is less 
than 2% which indicates the numerical uncertainties of the 
simulation results. 
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C. Static properties of bead-spring polymers 

In the next two sections we present results from simulations 
of bead-spring polymers of various chain length N. Here we 
look at the equilibrium values of the radius of gyration (R 2 g ) 
and the end-to-end (R 2 e ) distance. These are given by 



<*2> = 



2N 2 



£< 



(34) 



1 r i i i i r 
: "<§> 

q _ _ "-ft. 



o 



0.1 - 




•9. 



and 

{R 2 ee) = ((rN-r 1 f}. 
Their theoretically predicted scaling behavior is 



(R*) cx (R 2 ee ) <x(N-l) 



2;/ 



(35) 



(36) 



In a good solvent, perturbation analysisLLJ predicts an expo- 
nent of v ?s 0.588 for N — > oo, which is also reproduced 
in our simulations as shown in figure [lj For clarity, figure 
[I] only reproduces the results obtained with Fixman's Cheby- 
shev approximation and with our TEA-HI. Using the original 
Cholesky factorization of Ermak and McCammon or no HI at 
all gave results that were, within the numerical fluctuations, 
indistinguishable from the ones shown. However, for runtime 
reasons we only ran simulations for N < 70 with Ermak's 
original method. 

Over the range of chain lengths TV = 4 . . . 200, we ob- 
tained a ratio of {R^ e )/{R 2 gL~ 6.8, which is about 10% larger 
than the,*esults of Li et al.EZl, Jendrejack et al.t-J, or Liu and 
Diinwegt-J. The difference may be due to the fact that in our 
simulations the ratio between bond length L and bead radius 
a is smaller than used there. Consequently, the polymer can 
not be compacted as much as with relatively smaller beads. 

Even as the correct (R 2 ) and (R 2 e ) do not directly prove 
that our truncated expansion HI is correct, they show that it 
does not introduce any static perturbations to the polymer con- 
formations. 



D. Dynamic measures of bead-spring polymers 

A central dynamic property of a diffusing polymer is its 
center of mass diffusion coefficient D cm , which is predicted 
to scale as D cm cx N~ v with v w 0.588 with Rotne-Prager 
HI and oc 1/N without HI. The scaling without HI was re- 
produced in our simulations, as can be seen in figure ^|, which 
gives D cm in units of the diffusion coefficient Dq of an in- 
dividual bead. Without HI, the results are only shown up to 
N = 100, but were calculated for N = 2 . . . 400. 

With HI, the results with the original method of Ermak 
and McCammon and with Fixman's approximation for N = 
2 ... 70 were indistinguishable within the numerical uncer- 
tainties. For N > 70, only the faster Chebyshev approxi- 
mation was used. D cm from these simulation can be fitted 
well with A~ - 56 , which means that in our simulations the 
diffusion of long polymers was slightly faster than predicted 
theoretically. 
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FIG. 2: Diffusion coefficient D cm of the center of mass of bead- 
spring polymers of various chain lengths N from simulations without 
HI (diamonds) and with HI using Fixman's Chebyshev approxima- 
tion (open points) or our truncated expansion approximation (filled 
points). The dashed line indicates the observed scaling of jv _0 j6 . 
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FIG. 3: Relaxation times r corr of the autocorrelation function of the 
end-to-end vector R ee of polymers of various lengths TV from BD 
simulations without HI (diamonds) and with HI according to Fix- 
man's Chebyshev approximation (open points) and our truncated ex- 
pansion approximation (filled points), respectively. The relaxation 
times were obtained by fitting a stretched exponential to the autocor- 
relation function. The dashed line indicates the scaling cx TV 3 " with 
v = 0.588. 



The results with our TEA-HI show the same scaling behav- 
ior, while D cm is slightly smaller than with the correct factor- 
ization of D. The relative deviation between the results with 
our approach vs. the Chebyshev approximation is about 5%, 
i.e., our truncated expansion reproduces about 95% of the ef- 
fect of the hydrodynamic correlation on the overall diffusive 
motion of the polymer. 

Another measure, which is sensitive to the internal dynam- 
ics of the polymer, is the autocorrelation function (R ee (t) ■ 
R ee (0)) of the end-to-end vector R ee , which decays exponen- 
tially with a time constant r corr (the Zimm time). Figure ^| 
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shows that the fitted relaxation times are proportional to N 3 " 
as predicted theoretically. With our TEA-HI, the relaxation is 
again slightly slower by some 10%, which indicates that our 
approximation takes most of the hydrodynamic correlation be- 
tween the beads into account. Without HI, the relaxation of 
R ee takes about one order of magnitude longer at TV = 100 
and even more for longer chains. 



E. Runtime considerations 

Without HI, one needs to evaluate at each timestep the 
forces from the N — 1 springs connecting the TV beads and 
also the TV (TV — l)/2 repulsive two-body interactions. For 
each of the TV beads a random displacement has to be chosen 
and then the beads are moved according to the external and 
the random forces. Consequently, the observed runtime T per 
timestep can be fitted with T = t + ti TV + t 2 TV (TV - 1) 
as shown in figure Q For one million timesteps on a 2.8 GHz 
Pentium 4 CPU, we obtained to = 0.2 seconds, ti = 0.95 sec- 
onds, and ti = 0.038 seconds. To obtain the runtimes, we ran 
each of the simulations for several minutes until the runtime 
for one million timesteps could be calculated with sufficient 
statistical accuracy. 

With HI, also the hydrodynamic tensor has to be set up, for 
which TV (TV — l)/2 distances between the beads have to be 
calculated. For our truncated expansion, the TV normalization 
constants d ( |20| ) are evaluated in 0(N 2 ) time. For the ex- 
pansion coefficients fy, the average coupling e = (D^/Da) 
is required, for which the (TV — 1) (TV — 2)/2 off-diagonal en- 
tries of the diffusion tensor have to be summed up. Finally, 
the effective hydrodynamically corrected forces of equations 
( p"5| ) and ([16]) are evaluated in 0{N 2 ) time, leading to an over- 
all quadratic runtime that can be fitted with to = 1 seconds, 
ti = 0.2 seconds, and t 2 = 0.3 seconds for one million 
timesteps. For the very simple polymer model used, the in- 
clusion of HI thus slowed down the simulations by a constant 
factor of about ten for large TV. For more realistic — and thus 
more expensive — interactions between the individual beads, 
the relative cost of considering HI will even decrease. ._. 

In the original formulation of Ermak and McCammonU, 
the random displacements are determined via a Cholesky fac- 
torization of the diffusion tensor. This led to the expected 
runtime behavior of C(TV 3 ) in our simulations, which was fit- 
ted with T = 0.027TV 3 + 1.5TV 2 + 0.8TV seconds for one 
million timesteps (see figure Q). The more efficient Cheby- 
shev expansion led to an effective runtime of C(TV 2 5 ), i.e., 
T = 0.14TV 2 ' 5 4- 4TV at an accuracy of e f = 10" 3 . Reduc- 
ing ef to 10~ 2 or 10 _1 led to a small speedup of a factor 
of 2 or 3, respectively. We note that the Chebyshev expan- 
sion needs a lot of computer memory for the repeated matrix- 
vector multiplications, especially for large particle numbers. 
Consequently, at TV = 1000 the Chebyshev approximation 
was slowed down due to memory constraints of the computer 
we used and the runtime for TV = 2000 could not be deter- 
mined any more on this machine. For our approximation, the 
memory requirements were determined by the diffusion ten- 
sor, while the Cholesky factorization needed additional tem- 




10 100 1000 

N 



FIG. 4: Comparison of the runtime behavior of BD simulations with 
the different methods to include HI vs. the chain length TV of the 
bead-spring polymer. The datapoints give the times required to sim- 
ulate one million timesteps, the lines are polynomial fits to the data 
as explained in the text. 



porary storage of about the same size as the diffusion tensor. 

With the very simple interactions between the beads that we 
used here in the simulations, the inclusion of hydrodynamics 
with our approximation slowed down the simulation by a con- 
stant factor of about ten. A simulation of a polymer of 1000 
beads would consequently run for about half of a day with our 
approximated HI for every hour it takes without HI. The same 
simulation would take four to five days with Fixman's Cheby- 
shev approximation and more than a month with the original 
algorithm of Ermak and McCammon. 

The propagation was performed with the simple Euler prop- 
agator and a rather conservative timestep. The simulations 
therefore could be easily accelerated by using a more ad- 
vanced propagation scheme like, e.g., the semiimplicit scheme 
ofJendrejackEj or the Trotter expansion of De Fabritiis et 
alH 

As our approximation calculates the effective random 
forces from individual two-body contributions, it is now also 
possible to introduce distance dependent cut-offs to the hydro- 
dynamic interactions to reduce the number of terms that have 
to be summed up. Similar to the cut-offs for direct interac- 
tions, this will speed up the simulation. With the traditional 
methods of Ermak and McCammon or of Fixman, a cut-off 
would lead to vanishing entries in the diffusion tensor without 
reducing its size and, thus, the numerical effort to factorize 
it. With these algorithms, cut-offs only introduce errors in 
the dynamics without any gain. We did not consider cut-offs 
here, as it is a project on its own to investigate the trade-off be- 
tween the perturbations of the hydrodynamic interactions vs. 
the achieved runtime savings. 

Finally, we note that with the two-body contributions used 
in our TEA-HI, the hydrodynamic coupling can be summed 
up in parallel to the inter-particle forces. Then, there is no 
need to build up and keep the complete 3TV x 3TV hydrody- 
namic tensor in computer memory, only the actually required 
3x3 submatrices Dy of equations (^|) and (fy have to be de- 
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termined temporarily. Then the required memory for many- 
particle simulations increases only linearly with N. Conse- 
quently, our truncated expansion is not only faster than the 
originally proposed Cholesky factorization and the improved 
Chebyshev expansion, but it also needs much less computer 
memory, which allows for even larger systems to be handled 
efficiently. 

IV. CONCLUSIONS 

We showed how the correlated random displacements can 
be approximated efficiently in Brownian dynamics simula- 
tions with hydrodynamic interactions. In our truncated ex- 
pansion approximation, effective hydrodynamically corrected 
random forces are determined with the aim to reproduce the 
statistical moments of the thermal motion. Truncating the ex- 
pansion at the second order allows to calculate the random 
forces from individual two-body contributions in 0(N 2 ) time. 
For these only 3x3 submatrices of the 3iV x 37V hydro- 
dynamic tensor are required temporarily. Our approximation 
consequently has the same runtime and storage scaling as the 
calculation of the direct interactions between the particles. 

We then compared the simulation results with our approx- 
imation to the results obtained with the exact method of Er- 
mak and McCammon and with the Chebyshev approxima- 
tion of Fixman, the runtimes of which scaled as 0(N 3 ) and 



0(N 2 5 ), respectively. 

Both the dynamics of a bead-spring dimer with variable 
separation of the beads and of polymers of chain lengths of 
N = 2 . . . 200 showed that our approximation captures about 
95% of the hydrodynamically introduced correlation. For 
most applications on chemical or biological systems, where 
the direct interactions have to be approximated anyway, this 
appears a completely sufficient level of accuracy. Moreover, 
it allows for accounting of important hydrodynamic effects in 
systems where they were sofar mostly ignored for computa- 
tional reasons. 

As already mentioned above, we see the main application of 
our truncated expansion hydrodynamics in Brownian dynam- 
ics simulations of large biological systems, where one is in- 
terested in the dynamics of many-particle association.!, trans- 
port processes, or protein folding where the complicated in- 
teractions between the proteins or parts thereof have to be ap- 
proximated rather crudely. For these applications, it is surely 
better to include most of the effects of hydrodynamics with- 
out slowing down the simulation too much than to either have 
no HI at all or the correct HI at a prohibitive computational 
cost. Other applications could be to describe the behavior 
of non- | soherical particles by assembling them from smaller 
spheresOLil. Here the quality of the model increases with the 
number of spheres used to build the non-spherical particles. 
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